% TAB_OPT2.M
%  samples from the posterior density using Metropolis-Hastings
%  Proposal Density : Tailored Multiple-block Methods from Chib and Ramamurthy (2009)  
%
%
% Taeyoung Doh
% created: 03/01/2011
%--------------------------------------------------------------------
% house cleaning
%--------------------------------------------------------------------
clear all;
close all;
diary off;


% MHRUN

%	 1 for posterior draws
	
mhrun = 1;

%--------------------------------------------------------------------
% default parameters
%--------------------------------------------------------------------

append  = 1;		% append mode



% scaling factors, number of blocks 

    sbar   = 0.8; 
	nu     = 15;    % degree of freedom for t-distribution 
	nblock = 10;	% number of simulation blocks 
	nsimul = 100;	% number of simulation draws in each block


%--------------------------------------------------------------------
% model descriptions and basic calculations
%--------------------------------------------------------------------

% call environments and model selection
rsmodel_spec;
global msel 



lfile = [path_.para 'stat_mhrun21.mat']; 
load(lfile); 
  pmod = drawmean';
path_.mhrun = [path_.post 'mhrunpm_tab' num2str(mhrun)  '/']; 

chkdir(path_.mhrun);  
 
% global variables for parameter specification

% parameter names
pnames = feval(fnames_.pnames,msel);

%--------------------------------------------------------------------
% prior specification
%--------------------------------------------------------------------
global prior_;
prior_ = feval(fnames_.prior,msel);

%--------------------------------------------------------------------
% data and parameters
%--------------------------------------------------------------------

% load data
loaddata;
global data 
global paraspec_;
paraspec_ = feval(fnames_.prtr,msel);

% starting (or true) paremeters
para = pmod.*prior_.maskinv + prior_.fix.*prior_.mask;
nobs  = size(data,1);	% number of observations
npara = length(para);	% number of parameters


%--------------------------------------------------------------------
% starting point
%--------------------------------------------------------------------
if append		% append mode

	% load para_app.mat
	lfile = [path_.mhrun 'para_app.mat'];

	if exist(lfile,'file')
		load(lfile);

		paranew = paraapp';
		postnew = postapp;
		likenew = likeapp;

		if draws==nsimul
			iblock0 = blocks+1;
			iblockn = nblock;
			j0      = 1;
		elseif draws<nsimul
			iblock0 = blocks;
			iblockn = nblock;
			j0      = draws+1;

			lfile = [ path_.mhrun 'mhpm' num2str(iblock0,'%04.0f') '.mat'];
			load(lfile);
		end

	
		% print cotinuing point
		fmt='%15.8f ';
		disp(' ');
		disp(' ');
		disp(sprintf('Metropolis-Hastings continues at ... '));

	else
		append  = 0;
    end   
end

  if ~append   
    iblock0 = 1;
	iblockn = nblock;
	j0      = 1;


paranew = para;
[postnew,likenew]=rsfnmh(para,data); 
  end 
% print starting point	% print head materials
		fmt='%15.8f ';
		disp(' ');
		disp(' ');
	disp(sprintf('----------------------------------------------------------------------'));
disp('Parameter  Fixed          Values  ');
for i=1:length(para)
	disp(sprintf(['%s      %d '  repmat(fmt,1,1) ],pnames(i,:),prior_.mask(i),paranew(i)));
end
disp(' ');
disp(sprintf(['Log Prior      : ' fmt],postnew-likenew));
disp(sprintf(['Log Likelihood : ' fmt],likenew));
disp(sprintf(['Log Posterior  : ' fmt],postnew));
disp(sprintf('----------------------------------------------------------------------'));
disp(' ');


%--------------------------------------------------------------------
% Metropolis-Hastings simulation at a starting point, PARANEW
%--------------------------------------------------------------------

% set clock
t0 = now;

% update
likeold = likenew;
postold = postnew;
paraold = paranew;

% block loop
for iblock=iblock0:iblockn;
      
	% set time
	t1 = now;

	% initialize block output
	if j0==1
		likesim = zeros(nsimul,1);
		postsim = zeros(nsimul,1);
		parasim = zeros(nsimul,npara);
	end

	% within-block loop
  for j=j0:nsimul
     
        rindx = randperm(npara);
        sub_ind = 1;
        s = rand(npara,1); 
        while sub_ind<=npara
          
          stop=0; 
          sub_block=rindx(sub_ind); 
           while (stop==0)    
                    sub_block=sub_block;       
                    if (s(sub_ind)<=sbar) && (sub_ind<npara)  
                     sub_ind=sub_ind+1;
                     sub_block=[sub_block;rindx(sub_ind)]; 
                    elseif (s(sub_ind)>sbar) && (sub_ind<npara) 
                     stop=1;    
                     sub_ind=sub_ind+1;
                    elseif sub_ind>=npara 
                     stop=1;   
                     sub_ind=sub_ind+1; 
                    end  
           end
                
              global paraold sub_block
              c0 = paraold; 
              fun =@rsfnmhtab_simul;
              maxiter = 50*length(sub_block);   % 20*length(sub_block) up to 1000, 50*length(sub_block) from 1001 
              maxfun  = 100*length(sub_block);  % 40*length(sub_block) up to 1000, 100*length(sub_block) from 1001  
              c0 = invtrans(c0); 
              hybridopts = optimset('Display','iter','MaxFunEvals',maxfun,'MaxIter',maxiter);
              options = saoptimset('Display','iter','MaxFunEvals',maxfun,'MaxIter',maxiter);
              [c0,FVAL,EXITFLAG]=simulannealbnd(fun,c0,[],[],options);
              if EXITFLAG~=1 
              [c0,FVAL,EXITFLAG]=fminsearch(fun,c0,hybridopts);
              end 
              
              
              % find the maximum for each sub_block 
              
              c0=trans(c0);
              paranew=paraold;
              paranew(sub_block)=c0(sub_block); 
              [postnew,likenew]=rsfnmh(paranew,data); 
      
            if  postnew>postold     % accept
               disp('acceptance')
                paraold = paranew;
               postold = postnew;
               likeold = likenew; 
            end
 
               likesim(j)   = likeold;
			   postsim(j)   = postold;
			   parasim(j,:) = paraold';
            
           
         
      end % loop for sub_block
            
            
             modsize = 1;
	

			% save draws in output matrices
			 

		% print log
		if mod(j,modsize)==0  
			disp(' ');
			disp(' ');
			disp(' ');
			diary on;
			disp(' ');
			disp(sprintf('block: %2.0f / %2.0f        simulation steps: %5.0f / %5.0f',iblock,iblockn,j,nsimul));
			disp(sprintf('-----------------------------------------------------')); 
			disp(         'Parameters             : ');
			for i=1:length(para)
				disp(sprintf(['     %s           '  fmt ],pnames(i,:),parasim(j,i)));
			end
			disp(' ');
			disp(sprintf(['Likelihood             : ' fmt],sum(likesim)/j));
			disp(sprintf(['Posterior              : ' fmt],sum(postsim)/j));
			disp(' ');

			disp(        ['Elapsed Time (Block)   :  ' datestr(now-t1,13)]);
			disp(		 ['Elapsed Time (Total)   :  ' datestr(now-t0,'dd') ' days,  ' datestr(now-t0,13) ]);
			disp(sprintf('-----------------------------------------------------'));
			diary off;

			%save results in the meanwhile
			sname = [ path_.mhrun 'mhpm' num2str(iblock,'%04.0f') '.mat'];
			save(sname, 'likesim', 'postsim', 'parasim');

			% save for append mode
			blocks  = iblock;
			draws   = j;
			paraapp = parasim(j,:);
			likeapp = likesim(j);
			postapp = postsim(j);

			sfile = [path_.mhrun 'para_app.mat'];
			save(sfile,'paraapp','likeapp','postapp','blocks','draws');
		end

end	% within-block loop

	%save results for each block
	sname = [ path_.mhrun 'mhpm'  num2str(iblock,'%04.0f') '.mat'];
	save(sname, 'likesim', 'postsim','parasim');

	% save for append mode
	blocks  = iblock;
	draws   = nsimul;
	paraapp = parasim(end,:);
	likeapp = likesim(end);
	postapp = postsim(end);

	sfile = [path_.mhrun 'para_app.mat'];
	save(sfile,'paraapp','likeapp','postapp','blocks','draws');

	j0 = 1;
 end		% block loop



%--------------------------------------------------------------------
% save and print
%--------------------------------------------------------------------

% calculate elapsed time and print into log
%diary on;
disp(' ');
disp('-------------------------------------------------------');
disp(['Start       :  ' datestr(t0)]);
disp(['End         :  ' datestr(now)]);
disp(['Elapsed Time:  ' datestr(now-t0,'dd') ' days,  ' datestr(now-t0,13) ]);
disp('-------------------------------------------------------');
%diary off;

